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ABSTRACT 



The expansion and collision of two wind-blown superbubbles is investigated numerically. Our 
models go beyond previous simulations of molecular cloud formation from converging gas flows by 
exploring this process with realistic flow parameters, sizes and timescales. The superbubbles are 
blown by time-dependent winds and supernova explosions, calculated from population synthesis 
models. They expand into a uniform or turbulent diffuse medium. 

We find that dense, cold gas clumps and filaments form naturally in the compressed collision 
zone of the two superbubbles. Their shapes resemble the elongated, irregular structure of observed 
cold, molecular gas filaments and clumps. At the end of the simulations, between 65 and 80 
percent of the total gas mass in our simulation box is contained in these structures. 

The clumps are found in a variety of physical states, ranging from pressure equilibrium with 
the surrounding medium to highly under-pressured clumps with large irregular internal motions 
and structures which are rotationally supported. 



Subject headings: ISM: bubbles - 
structure — stars: formation 

1. Introduction 



ISM: clouds — ISM: general — ISM: kinematics and dynamics — ISM: 



Even though molecular gas amounts to only a 
small fraction of the volume of the Galactic in- 
terstellar medm m, it dominates interstellar mass 
(jFerrierd 120011 ). Moreover, being the site of all 
star formation, it plays an essential role in galaxy 
evolution. 

Molecular gas exists in galaxies in highly irreg- 
ular clouds, with very clumpy structure and large 
non-thermal line widths, an indication of inter- 



(ISchneider fc Elmegreer] ll979L lLada et al.l Il999t 
HartmaniJl2002l lLada et al.l 1200/1 iMolinari eHd 



Sj), 



a fact that can be attributed to their for- 



nal supersonic turbulence ([Falgarone & Phillips 
1990HWilliams et al.llioQoh . The large-scale struc- 
ture of these clouds is very filamentary, ac- 
cording to observations in different wavelengths 



mat ion process ( MverslEo09f ). 

A natural mechanism for assembling the large 
amounts of gas necessary for the formation 
of molecular clouds is from large-scale con- 
verging atomic flows, such as expanding shells 
(|Elmegreen fc Ladall977tlMcCrav fc Kafatosll 19871 ) 
colliding shells ( Nigra et al.l 12008 ). or flows gen- 
erated by l arge-s c ale gravitationa l inst a bilitie s 



(|Yang et al.l 120071 : iKim fc Ostrikerl 120021 l2006f ) 



Gas in these collision environments is subject to 
a number of fluid instabilities that can cool and 
condense it enough for it to become molecular 
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(|Heitsch et al.ll2005f ). This scenario can not only 



explain the rapid formation of molecular hydrogen 
from atomic hydrogen, but also accounts for the 
turbulent nature of molecular clouds, attri buting 



it to their formation process 



2006). 



In the context of local star- formin g regions, the 



converging flow idea was proposed bv lBallesteros-Paredes et al 
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( 1999h and lHartmann et "all ( 2001 ) as an explana- 
tion of the fact that stellar populations in many 
star-forming clouds have age spreads that are sig- 
nificantly smaller than the lateral crossing time of 
the cloud. 

(|2001f ) proposed an alternative 



Pringle et al. 



mechanism for molecular cloud formation by ag- 
glomeration of gas which is already molecular, ar- 
guing that, in order for the resulting clouds to be 
shielded from the background UV radiation which 
could evaporate them, the material in the converg- 
ing flows had to be already dense and cold enough 
to be mostly molecular. Although the forma- 
tion of large clouds from smaller, already molec- 
ular clumps is a possible scenario ([Dobbsl ([20081 ), 
Dobbs, Pringle & Burkert, submitted), there is 
still no observational evidence for the existence of 
a large molecular gas reservoir in the ISM outside 



giant molecular clouds. Moreover, iBergin et al. 
(|2004l ) showed that dust shielding could allow gas 
of much smaller densities to form molecules behind 



shocks and Heitsch fe Hartmannl ( 2008 ) showed in 



colliding flows simulations with gravity that the 
collapse along the second and third dimensions can 
rapidly lead to the high column densities required 
for molecule formation and for efficient shielding 
from UV radiation. 

Heitsch et all ([20061 . l2QQ8bh studied the forma- 



tion of cold and dense clumps between two infinite 
flows which collide on a perturbed interface. Their 
study showed that cold structure can arise even 
from initially uniform flows if the conditions favor 
certain fluid instabilities, s uch as the Non- linear 
Thin shell instability (NTSI,|VisimiacJ (|l994h ). the 
Thermal Instability ( Fieldl fl965[ ) and the Kelvin- 
Helmholtz instability. 

used col- 



Vazquez- Semadeni et al.l (|2 



liding cylindrical flows of tens of parsecs length, 
adding random velocity perturbations to the aver- 
age flow velocity. This allowed them to study star 
formation efficiencies and the fates of individual 
clouds. 

These numerical experiments have investigated 
the mechanisms which lead to the condensa- 
tion and cooling of atomic to molecular gas and 
which cause the complex structure of the result- 
ing clouds, independent of the specific mechanism 
driving the flows. In a further improvement of 
these models, we study the possibility for molecu- 
lar cloud formation from colliding flows of limited 



thickness, such as thin shells from expanding su- 
perbubbles, and without explicit control of the 
collision parameters. 

We present two-dimensional hydrodynamical 
simulations of colliding superbubbles, both in a 
uniform and in a turbulent diffuse medium. The 
superbubbles are produced by stellar feedback in 
the form of time-varying winds and supernova 
explosions, calculated from population synthesis 
models. We show that in this environment per- 
turbations in the colliding flows arise naturally 
through the growth of instabilities at the edge of 
the bubbles, forming cold and dense clumps with 
a very filamentary structure and high internal ve- 
locity dispersions. We also find that, although tur- 
bulence naturally arises even in an initially quies- 
cent background, the presence of background tur- 
bulence introduces anisotropics in the shell frag- 
mentation and more pronounced filaments. 

In section [2] we explain our numerical method, 
in section [3] we present the results of our simula- 
tions and in section |4] we discuss our results and 
their implications. 

2. Numerical Method 
2.1. Numerical Setup 

We perform two-dimensional hydrod ynamical 
simulations using the RAMSES code (iTevssierl 
20021 ). which uses a second-order Godunov scheme 
to solve the equations of ideal hydrodynamics on 
a Cartesian grid. In these calculations we have 
not made use of the Adaptive Mesh Refinement 
(AMR) feature of the code. 

We simulate a region of physical size equal to 
500 2 pc 2 , with a resolution of 4096 points at each 
dimension. Thus we achieve a spatial resolution 
of about 0.1 pc. We have also performed an 8 192 2 
resolution run but, due to its high computational 
cost, only for limited integrations and for test- 
ing some resolution effects. Although a 0.1 pc 
resolution is not sufficient to resolve the small- 
est struct ures in the ISM (see for exa mple, dis- 



Hennebelle fc Auditl (|2007ah ). we can 



cussion m 

resolve many of the clumps which form in these 
simulations adequately. In our analysis we do not 
take into account structures that fall near our res- 
olution limit (see discussion in section [3]). 

As an initial condition we use either a uniform 
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or a turbulent diffuse medium, of hydrogen density 
tih = 1 /cm 3 and temperature T=8000 K. In the 
case of the turbulent medium these are just aver- 
age quantities. We assume an ideal monoatomic 
gas with a ratio of specific heats 7 equal to 5/3 and 
mean molecular weight \i = 1.2m#. The setup of 
the turbulent medium is discussed in more detail 
in section [231 

In this diffuse medium we set up time-dependent 
winds that are meant to mimic the combined ef- 
fects of winds and supernova explosions in young 
OB asociations (- see section [2^21 for details). We 
include two wind regions, placed on the edges of 
the computational domain, that is, 500 pc apart, 
and we assume them to form simultaneously. We 
follow their expansion into the diffuse medium 
until they collide and a turbulent region arises at 
their interaction region. Reflecting boundaries are 
used along the x-direction and outflow boundaries 
along the y-direction. 

O u r runs include a cooli ng;- h eating; function fol- 
l owing I Wo Hire et al.l (|l995f ) and lDalgarno fc McCrav 
(|l972[ ) for t he low temperature regime (T < 
25000if ) and ISutherland fc Dopital (|l993l ) in the 
high-temperature (25000if < T < 10 8 K) regime. 
The metallicity assumed in the simulations is so- 
lar. Our initial condition is chosen to lie on a sta- 
ble point of the heating-cooling equilibrium curve, 
so that no cooling instability can occur unless the 
medium is externally perturbed. Of course, this is 
not exactly true everywhere for the turbulent case, 
but density fluctuations are in any case not large 
enough to lead to a two-phase medium without 
triggering from the bubbles. 

These calculations include no external gravity 
field or self-gravity of the gas. 

Time zero for the simulation is when star for- 
mation occurs in the wind areas. The system is 
evolved in time until boundary effects start be- 
coming potentially important, which is 7 Myrs for 
both simulations. 

2.2. Wind model 

We simulate the mass and energy input from a 
young OB association t o its surroundings accord- 
ing to Voss et all (|2009f ). The population synthe- 
sis model presented there provides the energy and 
mass injection from an average star of an OB as- 
sociation with time. 



Figure [T] shows the energy and mass injection 
rate from one such star with time, from the com- 
bined effects of stellar winds and supernova explo- 
sions. More details on the wind implementation in 
RAMSES can be found in Fierlinger et al. (2010, 
in prep.) 

Two circular regions are selected, located on the 
edges of our computational domain and 50 of these 
average stars are placed in each domain. These 
stars insert energy and momentum to the selected 
region at each timestep, according to the model 
illustrated in figure [TJ Although we would be able 
to resolve the spatial distribution of stars, we avoid 
this complication and treat the whole population 
as a uniform source, occupying a region of radius 
equal to 10 parsecs. This also minimizes the initial 
noise in the expanding shell. 

2.3. Turbulence setup 

In order to achieve a turbulent initial condi- 
tion for the diffuse mediu m, we set up a turbulen t 
velocity field according to Mac Low et al.l ([ 19981 ). 
that is, we introduce Gaussian random perturba- 
tions in Fourier space, in a range of wave numbers 
from k = 1 to k = 4. 

The rms Mach number of the turbulent flow 
is chosen close to unity, so that the turbulent 
kinetic energy of the gas is equal to its ther- 
mal energy. The energy equipartition assump- 
tion between thermal and turbulent kinetic energy 
is consistent with the turbulent velocity disper- 
sions calculated for Ga lactic HI ( Verschuur 2004t 
Haud fc Kalberlall2007l ). 

After this initial velocity field is calculated, it 
is introduced as an initial condition to the code, 
using a uniform hydrogen density of n# = 1/cm 3 
and temperature T=8000 K and follow its evo- 
lution isothermally and with periodic boundary 
conditions until the density-weighted power spec- 
trum of the turb ulent velocity field has a power- 
law slope close to lKolmogorovl ()194lh and the den- 
sity field loses the signature of the initial condi- 
tion. During this time the Mach number is kept 
constant by driving. 

The resulting velocity, density and pressure 
structure is used as an initial condition for the 
simulation of bubble expansion, which is otherwise 
identical to the simulations in a uniform medium. 

During the calculations we have neglected the 
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driving of turbulence in the diffuse medium, since 
the turbulence crossing time for our computational 
domain, t cr c± 86 Myrs is much longer than our en- 
tire calculation. This makes driving unnecessary 
during the simulation, as there are no significant 
energy losses due to dissipation on this time scale. 

3. Results 

3.1. Simulations in a uniform background 
medium 

In order to study the expansion of the su- 
perbubbles and the development of fluid insta- 
bilities relevant to it, we first simulated the two 
wind regions in a homogeneous diffuse background 
medium. Figure [2] shows two snapshots of the sim- 
ulation in logarithm of temperature and logarithm 
of density. 

During the expansion of the superbubbles we 
observe three effects from three different fluid in- 
stabilities. The acceleration of the shock leads to 
the NTSI, which focuses material on fluctuation 
peaks. These condensations are unstable to Ther- 
mal Instability, so they condense and cool further. 
The velocity shear caused by the NTSI at the same 
time also triggers the Kelvin-Helmholtz instability. 

As one can see in figure [2j the NTSI develops 
faster along the x and y axes, where we observe 
more pronounced "finger- like" structures, charac- 
teristic of th is instability. Th i s is c learly a reso- 
lution effect. IVishniac fe Rvul (|l989h showed that 
for an expanding decelerating shock there is a crit- 
ical overdensity with respect to the post-shock gas 
above which the shell is unstable. The authors es- 
timated this critical overdensity to be of the order 
25 for a wind-blown shock. Since in our simula- 
tions we do not resolve the smallest cooling length 
adequately, gas cannot be condensed into as small 
a volume as it should be according to its cooling 
rate. Along the x and y axes though, for a Carte- 
sian grid, the grid cells are closer together than 
along the diagonals, leading to an effective smaller 
distance over which the shock can compress the 
gas in one timestep. Thus along the x and y 
axes the fluid can be compressed slightly faster, 
then the critical overdensity can be achieved ear- 
lier and the instability arises earlier. As mentioned 
in previous sections, we have performed a test run 
with 8192 2 resolution for comparison, but have 
not been able to suppress this feature. However, 



the faster growth of the instability along these two 
lines does not affect the average clump formation 
time and the clump properties in any measurable 
way. 

As mentioned above, the resulting cold and 
dense structures are a result of the Thermal In- 
stability. However, our simulations do not in- 
clude thermal conduction and therefore cannot, 
by definition, fulfill the "Field criterion". Simu- 
lations not including thermal conduction may be 
susceptible to artificial phenomena near the reso- 
lution limit with no apparent convergence with in- 
creasing resolution, since thermal conduction has 
a stabilizing effect against t he Thermal Instability 
(jKovama fe Inutsukall2Qo"3 ). 

Th is stabilization is discussed in lBurkert fc Lin 
who studied the Thermal Instability both 
analytically and numerically for a generic, power- 
law cooling function. They find that density per- 
turbations below a certain wavelength are always 
damped by thermal conduction. This wavelength 
depends on the shape of the cooling function, and 
in our case would be roughly 4 times the Field 
length. For typical warm ISM values, that is, a 
thermal conductivity of about 10 4 ergs cm -1 K -1 
sec -1 , a temperature of T=8000 K and a hydro- 
gen number density of n=l cm -3 , and for a cooling 
rate of 10 -26 ergs sec -1 , the Field length is about 
0.03 pc and four times this length is about our 
resolution limit. This means that, even though 
our resolution is marginal for resolving the small- 
est unstable fluctuations, higher resolution with- 
out thermal conduction would only produce arti- 
ficially small structure. Of course, for lower tem- 
peratures and higher densities the Field length is 
much smaller, so unstable density fluctuations due 
to Thermal Instability within our clumps are not 
resolved. 

Figure [3] shows a close-up of the last snapshot 
of this run, in logarithm of density and pressure. 
It is easier to identify the clumps in this figure, so 
one clearly sees all of them are located in under- 
pressured regions with respect to the rest of the 
gas. 

We identify clumps by selecting locations with 
hydrogen number density greater than 50 cm -3 
and temperature smaller than 100 K and using 
a friends-of friends algorithm to link such adja- 
cent locations together. The group of cells is then 
identified as a single clump. The density and tem- 
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Fig. 2. — Superbubble collision snapshots in a uniform diffuse medium. Plotted on the top panels is the 
logarithm of the hydrogen number density in log(cm -3 ) and on the bottom panels the logarithm of the gas 
temperature in log(K). Left: 3 Myrs after star formation, right: 7 Myrs after star formation 
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Fig. 3. — Zoom-in of the last snapshot (7 Myrs after star formation) of the uniform diffuse medium run. 
Plotted on the left is the logarithm of the hydrogen number density in log (cm -3 ) and on the right the 
logarithm of the gas pressure in log(K cm -3 ). The axes coordinates are in parsecs. 
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perature threshold is clearly an approximation, to 
account for the fact that our simulations do not 
include molecule formation and our cooling func- 
tion assumes optically thin gas. For our analysis 
we only use structures which contain more than 16 
cells, as smaller structures do not contain sufficient 
information and are considered under-resolved. 

Figure |4] shows the gas phase diagram at the 
end of the simulation. The plot shows all the gas 
in the simulation, color-coded to its correspond- 
ing mass fraction in the simulated box. The solid 
line is the cooling-heating equilibrium curve for 
the warm and cold gas and the red crosses are 
average clump properties. As an indication for 
the gas temperatures, three isotherms have been 
overplotted. Most of the gas mass seems to be in 
the cold phase (- see also figure [5j) , but there is 
also about 14% of the gas mass in the warm, ther- 
mally unstable regime. This gas is pushed from 
the stable to the unstable regime by the momen- 
tum inserted by the wind. As explained in section 
3.31 some of this thermally unstable gas is located 
around the cold clumps, forming a warm corona. 
The clumps lie on the equilibrium curve, as does 
the coldest g as, as expected from clumps formed 
by the Thermal Instability. 

Figure [5] shows the mass fraction of the com- 
putational domain in each phase of the gas with 
time. We define cold gas to have temperatures 
T < lOOif , warm gas to lie in the temperature 
regime of 100K < T < 25000K and hot gas to 
have temperatures T > 25000if . As expected, we 
start with almost entirely warm gas in the box and 
as time goes by more and more cold gas is created. 
Since the hot gas is very dilute, it only amounts 
to approximately 1 percent of the mass through- 
out the simulation. In the end of the simulation 
we have approximately 85 per cent of the mass in 
cold clumps and 14 percent in the warm phase. 
The mass injected by the OB association amounts 
to less than 10 ~ 5 of the total mass in the domain 
throughout the simulation. 

Figure [6] shows the total number of identified 
clumps with time for each simulation. At the end 
of the uniform background simulation, more than 
400 clumps have formed. This is not only caused 
by the TI creating new clumps, but also by the 
fragmentation of already existing clumps. 
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Fig. 1. — Time dependence of wind properties for 
an average star. The solid line shows the energy 
injection rate in units of solar luminosity and the 
dashed line shows the mass injection r ate in solar 



masse s per Myr. The data are from IVoss et al. 
(l2QQ9h 
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Fig. 4. — Pressure versus hydrogen density of the 
fluid in the box for the uniform medium run. The 
plot is from the last snapshot, 7 Myrs after star 
formation. The gas has been binned and color- 
coded according to the mass fraction it represents. 
The solid line is the cooling-heating equilibrium 
curve for the warm and cold gas and the dashed 
lines show the locations of three isotherms. Red 
crosses are average clump quantities. 
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Fig. 7. — Super-bubble collision snapshots in a turbulent diffuse medium. Plotted on the top panels is the 
logarithm of the hydrogen number density in log (cm -3 ) and on the bottom panels the logarithm of the gas 
temperature in log(K). Left: 3 Myrs after star formation, right: 7 Myrs after star formation. 
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Fig. 8. — Zoom-in of the last snapshot (7 Myrs after star formation) of the turbulent diffuse medium run. 
Plotted on the left is the logarithm of the hydrogen number density in log (cm -3 ) and on the right the 
logarithm of the gas pressure in log(K cm -3 ). The axes coordinates are in parsecs. 
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3.2. Simulations with background turbu- 
lence 




2 4 6 

time (Myrs) 

Fig. 5. — Gas fractions with time for the uni- 
form medium run. The dashed curve shows the 
mass fraction in the warm gas phase, the solid 
curve shows the mass fraction in the hot gas phase 
and the dotted curve the mass fraction in the cold 
phase. (Hot phase: gas with T > 25000if , warm 
phase: gas with 100 K < T < 25000K , cold phase: 
gas with T < 100K) 
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Fig. 6. — Number of identified clumps with time. 
Clumps were only counted if they included more 
than 16 dense and cold cells. The solid curve cor- 
responds to the run in a uniform diffuse medium, 
while the dashed curve to the run in a turbulent 
diffuse medium. 



As for the uniform medium case presented in 
the previous section, here as well we observe the 
formation of cold and dense structures due to 
the combined action of the NTSI, the Kelvin- 
Helmholtz and the Thermal Instability. The main 
difference in this case is the anisotropy caused by 
the turbulent density and velocity field to the for- 
mation of clumps, as illustrated in figure [71 The 
shell on the left-hand side fragments very early on 
in some locations, already at less than 1 Myr af- 
ter star formation, but the shell on the right-hand 
side fragments much later, at around 3 Myrs after 
star formation. Rotating the initial conditions by 
180 degrees produces the exactly opposite effect. 

This is caused by the difference in the back- 
ground velocity field and is related to the range 
in wh ich the Vishniac instability exists (jVishniac 
19941 ). The smallest and largest unstable wave- 



length both depend on the relative velocity of the 
shock-bounded slab, in our case the thin expand- 
ing shell, and the background medium. Perform- 
ing a simple test with a uniform background veloc- 
ity field indicated that the instability indeed grows 
first where the relative velocity of the shell with 
respect to the background medium is largest, in 
agreement with Vishniac's analysis. 

Figure [8] shows a zoom-in of the last snapshot, 
showing a picture very similar to that in figure 
[3j where the different gas phases are almost in 
pressure equilibrium, with clumps located in re- 
gions with a pressure almost an order of magni- 
tude lower with respect to their surroundings. 

After the NTSI starts to grow, the condensed 
material immediately becomes thermally unstable 
and we see the formation of cold and dense clumps, 
as in the case of the uniform medium. However, 
due to the background velocity field, as explained 
above, there are now regions where the NTSI has 
a faster growth rate with respect to that of a static 
medium. This leads to a larger Kelvin-Helmholtz 
shear, which is the reason why in this simulation 
longer filaments are formed compared to the uni- 
form background medium simulation. Some of 
these structures are illustrated in figures [lOj [TT] 
and [I2j These figures are zoomed-in fractions of 
the full domain, shown in figure [9] 

Following these filaments from their formation 
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Fig. 12. — Time evolution of a single filament. The plots show logarithm of hydrogen number density. The 
black contour shows the level nn = 50 cm -3 . From top left to bottom right: 4.3, 4.6, 5.3 and 7 Myrs after 
star formation. Note the change in scale between the snapshots. 



11 



log n H (cm" 3 ) 
^ — : : : 1 : i 

-4.00 -2.72 -1.43 -0.15 1.13 2.42 3.70 








100 200 300 400 500 



Fig. 9. — Simulation snapshot at t=5.3 Myrs af- 
ter star formation. Plotted here is the logarithm 
of the hydrogen number density. The axes are 
marked in parsecs. The black rectangleds show 
the positions of the filaments shown in figures [TUJ 
[TTland the 5.3 Myrs snapshot of figure fT2l 
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Fig. 10. — A large filament containing several 
smaller clumps. Plotted here is the logarithm of 
the hydrogen number density in log(cm -3 ). The 
snapshot corresponds to the white box centered 
at x=250 pc, y=200 pc in the t=5.3 Myrs snap- 
shot shown in figure [9j The axes are marked in 
parsecs from the axes origin. The black contour 
corresponds to n# = 50cm -3 



to the end of the simulation, we observe that they 
increasingly become more elongated and they con- 
stantly fragment into smaller clumps. Figure fT2l 
shows the time evolution of one of the filaments. 
Before the shells collide, the structure is a small 
clump of cold gas. During the collision, this clump 
gets caught in a large-scale shear and becomes 
more and more elongated until, at the end of the 
simulation, it has reached a total length of about 
100 pc. Note that during the bubble expansion 
the clump is not located inside the hot bubble, 
but rather at the edge of the shock. 

The anisotropy in the growth rate of the NTSI 
leads to the formation of much less clumps in this 
simulation with respect to the uniform background 
simulation, as illustrated in figure [6] 

Figure [13] shows the phase diagram for the gas 
in this simulation, in the same way as figure |4] 
We note two main differences with respect to the 
uniform background medium run. First, that the 
pressure at low densities (n < 0.1 cm -3 ) is higher 
than in the uniform medium run. This effect can 
be attributed to the additional compression from 
the turbulence. Second, we note that the clumps 
in this simulation are at not only at lower pressures 
with respect to the hot medium than the clumps in 
the uniform background run, but are also at lower 
absolute pressures than those clumps. This can 
be explained by the formation mechanism of the 
clumps in these simulations. Structures created by 
the Thermal Instability have lower pressures than 
their surroundings. When Thermal Instability has 
stopped acting on them, they move towards pres- 
sure equilibrium with the surrounding medium in 
approximately a sound crossing time, c s /L, where 
c s the sound speed and L the typical size of these 
clumps. For a typical density of 10 ~ 21 g/cm 3 and 
a typical temperature of 50 K, the sound speed 
in these clumps is about 0.8 km/sec. At this ve- 
locity sound waves will cross a 2 pc long clump 
in about 0.6 Myrs. In the turbulent diffuse back- 
ground run, clumps are formed later than in the 
uniform background run, which means they are 
less evolved and farther from equilibrium. Accord- 
ing to the above calculations, the time difference 
between clump formation between the two simu- 
lations is about 3 sound crossing times, which is 
enough for many the clumps in the uniform diffuse 
background simulation to have reached approxi- 
mate pressure equilibrium with the hot medium. 
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Fig. 11. — As in figure fTUl The snapshot corre- 
sponds to the white box centered at x=225 pc, 
y=314 pc in the t=5.3 Myrs snapshot shown in 
figure [9l 
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Fig. 13. — Pressure versus hydrogen density of the 
fluid in the box for the turbulent medium run. The 
plot is from the last snapshot, 7 Myrs after star 
formation. The gas has been binned and color- 
coded according to the mass fraction it represents. 
The solid line is the cooling-heating equilibrium 
curve for the warm and cold gas and the dashed 
lines show the locations of three isotherms. Red 
crosses are average clump quantities. 



Figure [14] shows the evolution of the mass frac- 
tions of the gas in different temperature regimes. 
Unlike the uniform medium run, in this simulation 
there is a maximum in the hot gas mass fraction 
at about 3.5 Myrs after star formation and a later 
minimum of the warm gas fraction at around 5 
Myrs after star formation. This can be explained 
by the delayed formation of cold gas in this simu- 
lation. Until the cold gas is formed, the inserted 
hot gas just compresses the warm gas, so the hot 
gas mass fraction increases, while the warm gas 
mass fraction decreases. When cold gas starts to 
form, it quickly dominates the mass, causing the 
warm and hot mass fractions to drop. 

3.3. Clump properties 

As shown in figure [6j we form hundreds of 
clumps in each snapshot of the simulations. This 
makes it very difficult to study each of them in de- 
tail, especially since they are constantly merging 
and splitting. However, we can make some gen- 
eral comments on the morphology of the identified 
clumps and look into some of their properties. 

Clumps are generally at lower pressures with 
respect to their surrounding gas (- see figures fT5| 
[T6irm and IT8 |) . This means that Thermal Instabil- 
ity is still acting in these regions, causing clumps 
to condense further. Two examples of condensing 
clumps are shown in figure [T5l Others (about 11% 
in the last snapshot of the uniform background 
run and 16% in the turbulent background run) 
are rotating. Two examples of rotating structures 
are shown in figure HH Rotation is either com- 
bined with compression or it introduces a centrifu- 
gal force which makes the clump expand. Of the 
rotating cores in the last snapshots of the uniform 
and the turbulent background runs, about 25% 
and about 35% respectively are at the same time 
condensing and the rest are expanding. In very 
few cases, less than 1%, the centrifugal force ex- 
actly balances the force due to the pressure differ- 
ence between the interior and the exterior of the 
clump. 

A small fraction of the identified clumps are in 
pressure equilibrium with their surroundings, at 
least at their most central parts. These clumps 
do not tend to host significant internal motions (- 
Figure [[8]). 

All clumps are surrounded by a warm corona 



13 



which is more dilute and at an intermediate pres- 
sure between their pressure and the one of the sur- 
rounding gas. This corona usually surrounds more 
than one clump, indicating that they are parts of 
larger structures, such as the ones illustrated in 

figures HOI EH and CE2 

In the following we focus on the velocity disper- 
sions, sizes and possible evolution of the clumps. 

3.3.1. Velocity dispersions 

Figure [19] shows the internal velocity disper- 
sions of the clumps at different times for both 
simulations. The velocity dispersion here is de- 
fined as the square root of the variance among all 
locations which compose the clump. This means 
that for many clumps a velocity dispersion may 
indicate rotation, compression, expansion or ran- 
dom motions. As random motions here we define 
any combination of rotation and compression or 
expansion. Since the clumps have no significant 
internal density fluctuations, mass- weighted veloc- 
ity dispersions are not very different from the ones 
presented here. 

Although many of the formed clumps host su- 
personic motions in their interior at all times, 
there is an indication that these motions decrease 
at later times, as shown in the velocity disper- 
sion plot of the last snapshots. This effect seems 
slightly more pronounced for the turbulent back- 
ground run, where the lagest internal velocity dis- 
persions have disappeared in the last snapshot. 

In the case of compression or expansion, this 
means that the clumps gradually move to pres- 
sure equilibrium with their surrounding gas. In 
the case of random motions, it indicates that these 
motions are inherited by the turbulent environ- 
ment that created the clumps but, in absence of 
any mechanism to sustain them, they die out. In 
the case of rotation, though, the situation is a bit 
more complicated. Rotation could be an effect of 
the large-scale Kelvin-Helmholtz shear, it can orig- 
inate from Thermal Instability accretion or can be 
a result of structures splitting or merging. An- 
gular momentum conservation sustains these mo- 
tions for longer, so we would only expect them to 
decrease on a viscous timescale. 

In order to study if the decrease in internal 
velocity dispersion is observable for individual 
clumps, we tracked individual clumps back in pre- 
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Fig. 14. — Gas fractions with time for the tur- 
bulent medium run. The dashed curve shows the 
mass fraction in the warm gas phase, the solid 
curves shows the mass fraction in the hot gas phase 
and the dotted curve the mass fraction in the cold 
phase. (Hot phase: gas with T > 25000if , warm 
phase: gas with 100K < T < 25000K, cold phase: 
gas with T < 100K) 
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Fig. 20. — Time evolution of the velocity disper- 
sions over sound speed for three isolated cores. 
The green line corresponds to a condensing clump, 
the blue line to a rotating clump and the red line 
to a clump with random motions. The data for 
each clump begin from the snapshot where we can 
still identify the clump as being the same. 
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Fig. 15. — Two examples of condensing clouds. Top panels show the logarithm of hydrogen number density 
in log(cm -3 ) and bottom panels show the logarithm of thermal pressure in log(K cm -3 ). The black arrows 
show the velocity field with the mean velocity of the central clump subtracted. Overplotted in black are the 
contour levels for nn equal to 50 cm -3 , 100 era -3 and 1000 cm~ s . 
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Fig. 16. — Two examples of rotating clumps. Top panels show the logarithm of hydrogen number density 
in log(cm -3 ) and bottom panels show the logarithm of thermal pressure in log(K cm -3 ). The black arrows 
show the velocity field with the mean velocity of the central clump subtracted. Overplotted in black are the 
contour levels for nn equal to 50 cm -3 , 100 era -3 and 1000 cm~ s . 
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Fig. 17. — Two examples of clumps hosting random motions. Top panels show the logarithm of hydrogen 
number density in log(cm -3 ) and bottom panels show the logarithm of thermal pressure in log(K cm -3 ). 
The black arrows show the velocity field with the mean velocity of the central clump subtracted. Overplotted 
in black are the contour levels for uh equal to 50 cm -3 , 100 cm -3 and 1000 cm -3 . 
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Fig. 18. — Two examples of clumps with small internal velocities. Top panels show the logarithm of hydrogen 
number density in log(cm -3 ) and bottom panels show the logarithm of thermal pressure in log(K cm -3 ). 
The black arrows show the velocity field with the mean velocity of the central clump subtracted. Overplotted 
in black are the contour levels for uh equal to 50 cm -3 , 100 cm -3 and 1000 cm -3 . 
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Fig. 19. — Clump internal velocity dispersion versus their sound speed at different times. From left to right, 
3, 5 and 7 Myrs after star formation. The top panel corresponds to the uniform background medium run, 
the bottom panel to the turbulent background medium run. 



vious snapshots and plotted their velocity disper- 
sions with time. Since dense structures are cre- 
ated all the time and clumps merge or split at 
each snapshot, it is very difficult to construct an 
algorithm able to automatically identify the same 
clump in different snapshots. Instead we identified 
the clumps by eye according to their positions and 
translational velocities. We focus here on three 
examples. We selected a rotating clump, a con- 
tracting clump and a clump with random internal 
motions. Their velocity dispersions over their cor- 
responding sound speed as a function of time are 
shown in figure [20] Since clumps evolve almost 
isothermally, their sound speed does not vary sig- 
nificantly during the interval shown in the plot. 
The green line corresponds to a contracting clump, 
the blue line to a rotating clump and the red 
line to a clump with random motions. All three 
clumps show some decrease in velocity dispersion 
with time. Although we are not able to make a 
complete study at this stage, we observe that the 
rotating clump practically maintains the same ve- 
locity dispersion throughout its existence, showing 
only a slight decrease, while the condensing clump 
and the clump with random motions show a de- 
crease in velocity dispersion. For the condensing 
clump this decrease is very pronounced. 



3.3.2. Sizes 

Figure [21] shows the size distribution of the 
clumps at different times for both simulations. 
The sizes are calculated as the square root of the 
area occupied by the clump. As indicated by the 
figure, we do not form clumps larger than ap- 
proximately 3 pc. However, the maximum clump 
length can reach about 10 pc. This means that, 
although structures may have one very large di- 
mension, they occupy a very small area. 

Figure [22] shows the distributions of the ra- 
tios of clump sizes over their corresponding Jeans 
length, at different times. The size distributions 
of the clumps have a very similar shape and range 
between the two simulations. As clumps are fairly 
uniform in density and temperature, their Jeans 
length does not vary significantly within each sin- 
gle clump. Clumps in general seem to be smaller 
than their corresponding Jeans length but, as time 
advances, some clumps become large enough to be 
potentially unstable to gravitational collapse. Of 
course, as we have not included gravity in these 
simulations we cannot know if this would actually 
be the case. 

As mentioned earlier, since the formed struc- 
tures are very filamentary, they are likely to con- 
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Fig. 21. — Clump size distributions in parsecs at different times. The plot on the left-hand side corresponds 
to the run in a uniform diffuse medium and the plot on the right to the run in a turbulent diffuse medium. 
The solid black, dashed red and dash-dotted blue histograms correspond to the 3 Myrs, 5 Myrs and 7 Myrs, 
respectively. 
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Fig. 22. — Clump size over clump Jeans length distributions at different times. The plot on the left-hand 
side corresponds to the run in a uniform diffuse medium and the plot on the right to the run in a turbulent 
diffuse medium. The solid black, dashed red and dash-dotted blue histograms correspond to the 3 Myrs, 5 
Myrs and 7 Myrs, respectively. 
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tain more Jeans lengths along a single dimen- 
sion. Note also that the algorithm we use to find 
clumps favors the identification of the smallest 
possible structures as separate entities. As men- 
tioned above however, the clumps we identify are 
usually parts of larger structures which are dy- 
namically interacting or surrounded by a common 
warm and more diffuse corona. 

3.3.3. Clump evolution 

Although most of the clumps are in low- 
pressure regions in the simulation, there are some 
clumps with approximately the same pressure as 
their surrounding gas. All of the clumps are 
surrounded by an intermediate pressure corona, 
which is also thermally unstable. 

This, in combination with the fact that clumps 
show a tendency of decreasing their internal mo- 
tions with time and that clouds in pressure equi- 
librium tend to host smaller internal motions leads 
us to believe that there might be an evolution from 
clumps out of equilibrium, with strong internal 
motions to more quiescent clumps. 

Figure [23] shows the tracks of the same three 
clumps we traced back in time on the pressure- 
density diagram. The dotted line is the cooling- 
heating equilibrium curve. The clumps form in 
an area where cooling dominates, possibly from 
thermally unstable gas, at the left of the figure. As 
time goes by they move on the equilibrium curve 
and gradually increase their density and pressure, 
while staying on the curve. The timescale of this 
evolution is about 2-3 Myrs, which corresponds to 
about 4-5 clump sound crossing times. 

4. Summary and Discussion 

In this work we have studied the combined ef- 
fects of the Vishniac, the Kelvin-Helmholtz and 
the Thermal Instability on the expansion and col- 
lision of two superbubbles, created by stellar winds 
and supernova explosions from young OB associ- 
ations. 

We have presented high-resolution, two-dimensional 
simulations to show that, even without simulating 
the effects of gravity, the shells formed by stellar 
feedback will condense and fragment to form cold 
and dense filaments. We have modeled the shell 
collision in a static, uniform diffuse medium and 
in a turbulent diffuse medium and found that 



background turbulence not only introduces an 
anisotropy in the shell fragmentation, but also 
helps to produce more filamentary structures. 
Our results provide a pictu re of the ISM similar 



to 



Audit Hennebelld ([20051 ) and lHennebelle fc Audit 



(|2007bh . where the ISM phases are tightly in- 



terwoven, with sharp thermal interfaces between 
them. Of course, a more detailed comparison is 
not possible, since that work was done with a much 
higher resolution, representing a smaller region of 
the ISM and not including the hot phase. In our 
simulations, elongated cold structures sit in warm, 
thermally unstable coronas, submerged in a hot di- 
lute medium that the stellar winds and supernovae 
from young OB stars create. Clumps are dynami- 
cal entities, constantly merging and splitting and 
in general not in pressure equilibrium with their 
surrounding gas. In that sense, the picture we see 
in our simulations is very diffe rent from the clas- 
sical I SM model proposed by iMcKee fc Ostriker 



([19771 ). where the clouds are treated as quasi-static 
spheres. The structures we find are rather long 
thin filaments, very similar to the large "bl obby 
sheets" proposed bv lHeiles fc Trolandl ([20031 ). 

We identify individual clumps in the simula- 
tions by setting density and temperature thresh- 
olds and connecting adjacent locations in the sim- 
ulation domain which exceed these thresholds. 
This method favors the identification of the small- 
est cold and dense structures. The clumps we 
identify tend to be connected in groups within a 
common warm corona, forming long filaments and 
are mostly located in areas of lower pressure with 
respect to their surroundings. The latter causes 
many of them to condense to higher densities and 
smaller volumes, almost isothermally. Apart from 
condensing clumps we find rotating clumps and 
clumps hosting irregular motions. All these inter- 
nal motions are reflected in the velocity dispersion 
as supersonic internal motions. Approximately 11 
to 16% of the identified clumps are rotating and 
are interesting candidates for forming protostel- 
lar disks. However, rotation should be studied in 
three dimensions for more reliable results. 

Internal clump motions tend to decrease with 
time, as the clumps come closer to thermal pres- 
sure equilibrium with their warm surroundings. 
However, if gravity and star formation was in- 
cluded in our simulations, we would expect these 
motions to be sustained for longer times due to 
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gravitational collapse and feedback from the newly 
formed stars in their interior. 

Cold clumps in our simulations form by Ther- 
mal Instability condensation of the ambient dif- 
fuse medium when it is perturbed by the expand- 
ing shells. This, in combination with the fact 
that the mass injected by the OB associations is 
a very small fraction of the total mass in our do- 
main throughout the simulation, shows that the 
material forming the new clumps could not have 
been enriched directly by the supernova explo- 
sions. The material ejected from the OB associa- 
tions will probably enrich the interstellar matter 
on much larger timescales, after the hot bubbles 
have mixed with the diffuse medium. 

Thermally unstable gas amounts to about 8 
to 10% of the total gas mass in the last snap- 
shots of our simulations, far from the almost 
50% that is commonly detected in observations 
([Heiles fc Trolandl l2QQ3h or found in simulations 



of turbulen t, thermally bistable flows (jGazol et al 



200ll l2QQ5h . This is because our simulations are 



dominated by what we identify as CNM, that is, 
gas with temperatures lower than 100K. Due to 
our very high resolution, this gas is also very dense, 
reaching hydrogen number densities of the order of 
10 5 cm -3 . This gas would be mostly molecular, so 
it would not be identified as cold HI in observa- 
tions. In simulations with gravity we would not 
expect to encounter this issue, since gravity would 
eventually dominate the dynamics of the formed 
clumps and turn them into stars once they be- 
came dense enough, so thermal instability would 
no longer be responsible for their evolution at 
these densities. 

The simulations presented in this paper are 
only a first step to modeling triggered molecular 
cloud formation using physically motivated collid- 
ing flow parameters. We have not attempted a 
parameter study in this paper, but it will be the 
object of future work to study the effect of varying 
the distance between the superbubbles, the num- 
ber of OB stars creating the superbubbles and the 
metallicity of the gas. 

In this work we have also not taken into ac- 
count galactic shear or density stratification. As- 
suming a galactic rotation rate of 26 km sec _1 
kpc _1 , the relative shear in our computational 
box would be 13 km sec -1 . This velocity would 
have a crossing time of approximately 37 Myrs, 



which is much longer than what our simulations 
last. Density stratification with a scale height of 
H=150pc would cause the superbubbles to expand 
faster along the vertical direction, but since we 
are interested mostly in what happens at the su- 
perbubble collision interface, we can neglect this 
effect as a first approximation. 

Our calculations also ignore magnetic fields and 
gravity. We would expect magnetic fields to play 
a significant role in the dynamics of the problem, 
both during the expansion of the superbubbles and 
also during the more complex collision phase, if 
they had a preferred orientation, but in two dimen- 
sions we would not be able to model all relevant 
phenomena properly anyway. Gravity is essential 
for studying the evolution of the clumps and for es- 
timating their star forming efficiency. Combining 
with the modeling of a third dimension, it would 
give us useful mass estimates of the formed struc- 
tures. Moving to three dimensions and including 
gravity is work in process and will be presented in 
a future paper. 

We find the main limitations of our work to 
be the lack of resolution and bidimensionality. 
As we have pointed out, our simulations do not 
reach the resolution required to capture all the 
relevant physical processes on the smallest scales. 
We have accounted for small-scale effects by only 
studying clumps which contain more than 16 grid 
cells. Higher numerical resolution will certainly 
provide more insight on the number of formed 
clumps and their internal structure. Although 
numerical effects do introduce a thermal conduc- 
tion effect, explicitly modeling thermal conduc- 
tion would help set the minimum scale for formed 
clumps and achieve convergence with increasing 
resolution. On the other hand, if gravity is sim- 
ulated it will probably already become important 
at length scales larger than the Field length. 

The restriction of the presented models to 
two dimensions is emphasizing the formation 
of filaments. In three dimensions, these struc- 
tures could have a sheet-like morphology. Grav- 
ity could then be re sponsible for turning these 
sheets into filaments ( Burkert fc Hartmarml 2004t 



Hartmann fc Burkertl [20071 : iHeitsch et al.ll2008aF . 
Although combining an increase in resolution with 
studying the problem in three dimensions is nu- 
merically very expensive, it is probably feasible 
with the use of AMR. 
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Fig. 23. — Isolated core tracks on a pressure- 
density diagram. The green triangles correspond 
to a condensing clump, the blue crosses to a ro- 
tating clump and the red diamonds to a clump 
with random motions. The solid black line is part 
the cooling-heating equilibrium curve which was 
shown in figures 2] and [13j 
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